#include "ManagedLasoProvider.h"

/* Table of constant values */

static integer c__0 = 0;
static integer c__1 = 1;


/* *********************************************************************** */

 /* Subroutine */int SmartMathLibrary::LasoPack::ManagedLasoProvider::dlabfc_
   (integer *n, integer *nband, doublereal *a, doublereal *sigma, integer
   *number, integer *lde, doublereal *eigvec, integer *numl, integer *ldad,
   doublereal *atemp, doublereal *d__, doublereal *atol)
{
  /* System generated locals */
  integer a_dim1, a_offset, eigvec_dim1, eigvec_offset, atemp_dim1,
    atemp_offset, d_dim1, d_offset, i__1, i__2, i__3;
  doublereal d__1, d__2;

  /* Builtin functions */
  //double d_sign(doublereal *, doublereal *);

  /* Local variables */
  static doublereal zero[1];
  static integer i__, j, k, l, m;
  //extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
  // doublereal *, integer *), dswap_(integer *, doublereal *, integer 
  // *, doublereal *, integer *), daxpy_(integer *, doublereal *, 
  // doublereal *, integer *, doublereal *, integer *);
  static integer la, ld, kk, nb1, lpm;


  /*  THIS SUBROUTINE FACTORS (A-SIGMA*I) WHERE A IS A GIVEN BAND */
  /*  MATRIX AND SIGMA IS AN INPUT PARAMETER.  IT ALSO SOLVES ZERO */
  /*  OR MORE SYSTEMS OF LINEAR EQUATIONS.  IT RETURNS THE NUMBER */
  /*  OF EIGENVALUES OF A LESS THAN SIGMA BY COUNTING THE STURM */
  /*  SEQUENCE DURING THE FACTORIZATION.  TO OBTAIN THE STURM */
  /*  SEQUENCE COUNT WHILE ALLOWING NON-SYMMETRIC PIVOTING FOR */
  /*  STABILITY, THE CODE USES A GUPTA'S MULTIPLE PIVOTING */
  /*  ALGORITHM. */

  /*  FORMAL PARAMETERS */


  /*  LOCAL VARIABLES */


  /*  FUNCTIONS CALLED */


  /*  SUBROUTINES CALLED */

  /*     DAXPY, DCOPY, DSWAP */


  /*  INITIALIZE */

  /* Parameter adjustments */
  a_dim1 =  *nband;
  a_offset = 1+a_dim1 * 1;
  a -= a_offset;
  eigvec_dim1 =  *lde;
  eigvec_offset = 1+eigvec_dim1 * 1;
  eigvec -= eigvec_offset;
  d_dim1 =  *ldad;
  d_offset = 1+d_dim1 * 1;
  d__ -= d_offset;
  atemp_dim1 =  *ldad;
  atemp_offset = 1+atemp_dim1 * 1;
  atemp -= atemp_offset;

  /* Function Body */
  zero[0] = 0.;
  nb1 =  *nband - 1;
  *numl = 0;
  i__1 =  *ldad **nband;
  dcopy_(&i__1, zero, &c__0, &d__[d_offset], &c__1);

  /*   LOOP OVER COLUMNS OF A */

  i__1 =  *n;
  for (k = 1; k <= i__1; ++k)
  {

    /*   ADD A COLUMN OF A TO D */

    d__[ *nband +  *nband * d_dim1] = a[k *a_dim1 + 1] -  *sigma;
    m = min(k,  *nband) - 1;
    if (m == 0)
    {
      goto L20;
    }
    i__2 = m;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      la = k - i__;
      ld =  *nband - i__;
      d__[ld +  *nband * d_dim1] = a[i__ + 1+la * a_dim1];
      /* L10: */
    }

    L20: 
    /* Computing MIN */
    i__2 =  *n - k;
    m = min(i__2, nb1);
    if (m == 0)
    {
      goto L40;
    }
    i__2 = m;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      ld =  *nband + i__;
      d__[ld +  *nband * d_dim1] = a[i__ + 1+k * a_dim1];
      /* L30: */
    }

    /*   TERMINATE */

    L40: lpm = 1;
    if (nb1 == 0)
    {
      goto L70;
    }
    i__2 = nb1;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      l = k -  *nband + i__;
      if (d__[i__ +  *nband * d_dim1] == 0.)
      {
        goto L60;
      }
      if ((d__1 = d__[i__ + i__ * d_dim1], abs(d__1)) >= (d__2 = d__[i__ + 
        *nband * d_dim1], abs(d__2)))
      {
        goto L50;
      }
      if (d__[i__ +  *nband * d_dim1] < 0. && d__[i__ + i__ * d_dim1] < 0. ||
        d__[i__ +  *nband * d_dim1] > 0. && d__[i__ + i__ * d_dim1] >= 0.)
      {
        lpm =  - lpm;
      }
      i__3 =  *ldad - i__ + 1;
      dswap_(&i__3, &d__[i__ + i__ * d_dim1], &c__1, &d__[i__ +  *nband *
        d_dim1], &c__1);
      dswap_(number, &eigvec[l + eigvec_dim1], lde, &eigvec[k + eigvec_dim1],
        lde);
      L50: i__3 =  *ldad - i__;
      d__1 =  - d__[i__ +  *nband * d_dim1] / d__[i__ + i__ * d_dim1];
      daxpy_(&i__3, &d__1, &d__[i__ + 1+i__ * d_dim1], &c__1, &d__[i__ + 1+
        *nband * d_dim1], &c__1);
      d__1 =  - d__[i__ +  *nband * d_dim1] / d__[i__ + i__ * d_dim1];
      daxpy_(number, &d__1, &eigvec[l + eigvec_dim1], lde, &eigvec[k +
        eigvec_dim1], lde);
      L60: ;
    }

    /*  UPDATE STURM SEQUENCE COUNT */

    L70: if (d__[ *nband +  *nband * d_dim1] < 0.)
    {
      lpm =  - lpm;
    }
    if (lpm < 0)
    {
      ++(*numl);
    }
    if (k ==  *n)
    {
      goto L110;
    }

    /*   COPY FIRST COLUMN OF D INTO ATEMP */
    if (k <  *nband)
    {
      goto L80;
    }
    l = k - nb1;
    dcopy_(ldad, &d__[d_offset], &c__1, &atemp[l *atemp_dim1 + 1], &c__1);

    /*   SHIFT THE COLUMNS OF D OVER AND UP */

    if (nb1 == 0)
    {
      goto L100;
    }
    L80: i__2 = nb1;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      i__3 =  *ldad - i__;
      dcopy_(&i__3, &d__[i__ + 1+(i__ + 1) *d_dim1], &c__1, &d__[i__ + i__ *
        d_dim1], &c__1);
      d__[ *ldad + i__ * d_dim1] = 0.;
      /* L90: */
    }
    L100: ;
  }

  /*  TRANSFER D TO ATEMP */

  L110: i__1 =  *nband;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    l =  *n -  *nband + i__;
    i__2 =  *nband - i__ + 1;
    dcopy_(&i__2, &d__[i__ + i__ * d_dim1], &c__1, &atemp[l *atemp_dim1 + 1],
      &c__1);
    /* L120: */
  }

  /*   BACK SUBSTITUTION */

  if (*number == 0)
  {
    return 0;
  }
  i__1 =  *n;
  for (kk = 1; kk <= i__1; ++kk)
  {
    k =  *n - kk + 1;
    if ((d__1 = atemp[k *atemp_dim1 + 1], abs(d__1)) <=  *atol)
    {
      atemp[k *atemp_dim1 + 1] = d_sign(atol, &atemp[k *atemp_dim1 + 1]);
    }

    /* L130: */
    i__2 =  *number;
    for (i__ = 1; i__ <= i__2; ++i__)
    {
      eigvec[k + i__ * eigvec_dim1] /= atemp[k *atemp_dim1 + 1];
      m = min(*ldad, k) - 1;
      if (m == 0)
      {
        goto L150;
      }
      i__3 = m;
      for (j = 1; j <= i__3; ++j)
      {
        l = k - j;
        eigvec[l + i__ * eigvec_dim1] -= atemp[j + 1+l * atemp_dim1]*eigvec[k +
          i__ * eigvec_dim1];
        /* L140: */
      }
      L150: ;
    }
    /* L160: */
  }
  return 0;
} /* dlabfc_ */
